{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Example: Virtual circuit calculations\n", "\n", "---\n", "\n", "This example notebook will demonstrate how to construct **virtual circuits** (VCs) for a **MAST-U-like tokamak** equilibrium.\n", "\n", "### What are Virtual Circuits and why are they useful?\n", "\n", "VCs can be used to identify which (active) poloidal field (PF) coils have the most significant impact on a set of specified plasma shape parameters (which we will refer to henceforth as **targets**). The targets are equilibrium-related quantities such as the inner or outer midplane radii $R_{in}$ or $R_{out}$ (more will be listed later on). \n", "\n", "More formally, the **virtual circuit** (VC) matrix $V$ for a given equilibrium (and chosen set of shape targets and PF coils) is defined as\n", "\n", "$$ V = (S^T S)^{-1} S^T, $$\n", "\n", "where $S$ is the **shape** (Jacobian) matrix:\n", "\n", "$$ S_{i,j} = \\frac{\\partial T_i}{\\partial I_j}. $$\n", "\n", "Here, $T_i$ is the $i$ th target and $I_j$ is the current in the $j$ th PF coil. We note that $V$ is simply the Moore-Penrose pseudo-inverse of $S$. In FreeGSNKE, $S$ will be calculated using finite difference (see below). \n", "\n", "### How can these VCs be used?\n", "\n", "Once we know the VC matrix $V$ for a given equilibrium (and its associated target values $\\vec{T}$), we can specify a perturbation in the targets $\\vec{\\Delta T}$ and calculate the change in coil currents required to acheive the new targets. The shifts in coil currents can be found via:\n", "\n", "$$ \\vec{\\Delta I} = V \\vec{\\Delta T}. $$\n", "\n", "Using $\\vec{\\Delta I}$ we can perturb the coil currents in our equilibrium, re-solve the static forward Grad-Shafranov (GS) problem, and observe how the targets (call them $\\vec{T}_{new}$) have changed \n", "in the new equilibrium vs. the old targets $\\vec{T}$.\n", "\n", " ---\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Generate a starting equilibrums\n", "\n", "Firstly, we need an equilbirium to test the VCs on." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "# build machine\n", "from freegsnke import build_machine\n", "tokamak = build_machine.tokamak(\n", " active_coils_path=f\"../machine_configs/MAST-U/MAST-U_like_active_coils.pickle\",\n", " passive_coils_path=f\"../machine_configs/MAST-U/MAST-U_like_passive_coils.pickle\",\n", " limiter_path=f\"../machine_configs/MAST-U/MAST-U_like_limiter.pickle\",\n", " wall_path=f\"../machine_configs/MAST-U/MAST-U_like_wall.pickle\",\n", " magnetic_probe_path=f\"../machine_configs/MAST-U/MAST-U_like_magnetic_probes.pickle\",\n", ")\n", "\n", "# initialise the equilibrium\n", "from freegsnke import equilibrium_update\n", "eq = equilibrium_update.Equilibrium(\n", " tokamak=tokamak,\n", " Rmin=0.1, Rmax=2.0, # Radial range\n", " Zmin=-2.2, Zmax=2.2, # Vertical range\n", " nx=65, # Number of grid points in the radial direction (needs to be of the form (2**n + 1) with n being an integer)\n", " ny=129, # Number of grid points in the vertical direction (needs to be of the form (2**n + 1) with n being an integer)\n", " # psi=plasma_psi\n", ") \n", "\n", "# initialise the profiles\n", "from freegsnke.jtor_update import ConstrainPaxisIp\n", "profiles = ConstrainPaxisIp(\n", " eq=eq, # equilibrium object\n", " paxis=8e3, # profile object\n", " Ip=6e5, # plasma current\n", " fvac=0.5, # fvac = rB_{tor}\n", " alpha_m=1.8, # profile function parameter\n", " alpha_n=1.2 # profile function parameter\n", ")\n", "\n", "# load the nonlinear solver\n", "from freegsnke import GSstaticsolver\n", "GSStaticSolver = GSstaticsolver.NKGSsolver(eq) \n", "\n", "# set the coil currents\n", "import pickle\n", "with open('simple_diverted_currents_PaxisIp.pk', 'rb') as f:\n", " current_values = pickle.load(f)\n", "for key in current_values.keys():\n", " eq.tokamak.set_coil_current(key, current_values[key])\n", "\n", "# carry out the foward solve to find the equilibrium\n", "GSStaticSolver.solve(eq=eq, \n", " profiles=profiles, \n", " constrain=None, \n", " target_relative_tolerance=1e-4)\n", "\n", "# updates the plasma_psi (for later on)\n", "eq._updatePlasmaPsi(eq.plasma_psi)\n", "\n", "# plot the resulting equilbria \n", "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)\n", "ax1.grid(True, which='both')\n", "eq.plot(axis=ax1, show=False)\n", "eq.tokamak.plot(axis=ax1, show=False)\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Initialise the VC class\n", "\n", "We initialise the VC class and tell it to use the static solver we already initialised above. This will be used to repeatedly (and rapidly) to solve the static GS problem when calculating the finite differences for the shape matrix." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from freegsnke import virtual_circuits\n", "\n", "VCs = virtual_circuits.VirtualCircuitHandling()\n", "VCs.define_solver(GSStaticSolver)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Define shape targets\n", "\n", "Next we need to define the shape targets (i.e. our quantities of interest) that we wish to monitor. There are a number of shape target pre-defined in FreeGSNKE:\n", "- \"R_in\": inner midplane radius. \n", "- \"R_out\": outer midplane radius. \n", "- \"Rx_lower\": lower X-point (radial) position.\n", "- \"Zx_lower\": lower X-point (vertical) position.\n", "- \"Rx_upper\": upper X-point (radial) position.\n", "- \"Zx_upper\": upper X-point (vertical) position.\n", "- \"Rs_lower_outer\": lower strikepoint (radial) position.\n", "- \"Rs_upper_outer\": upper strikepoint (radial) position.\n", "\n", "Note that the following targets require additional options:\n", "- \"Rx_lower\": approx. radial position of the lower X-point (R,Z).\n", "- \"Zx_lower\": approx. vertical position of the lower X-point (R,Z).\n", "- \"Rx_upper\": approx. radial position of the upper X-point (R,Z).\n", "- \"Zx_upper\": approx. vertical position of the upper X-point (R,Z).\n", "- \"Rs_lower_outer\": approx. (R,Z) position of the lower outer strikepoint.\n", "- \"Rs_upper_outer\": approx. (R,Z) position of the upper outer strikepoint.\n", " \n", "We'll see how these options are defined in a moment. While they are not all strictly required, having them is advisable as the default calculations may return spurious values in some rare cases. \n", "\n", "More can be added (via a feature request), though we should say that it would need to be generic enough such that its definition is well-defined across different tokamak geometries and plasmas.\n", "\n", "There is the option to specify **custom** shape targets if these do not work for you---we will see this shortly." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# define the targets of interest and use the VC object to calculate their values for the equilibrium above\n", "targets = ['R_in', 'R_out', 'Rx_lower', 'Zx_lower', 'Rx_upper', 'Zx_upper', 'Rs_lower_outer', 'Rs_upper_outer']\n", "\n", "# define the target options in a dictionary (approx. (R,Z) locations of the X-points)\n", "# this helps identify the correct X-point in the code\n", "targets_options = dict()\n", "targets_options['Rx_lower'] = np.array([0.6, -1.1])\n", "targets_options['Zx_lower'] = np.array([0.6, -1.1])\n", "targets_options['Rx_upper'] = np.array([0.6, 1.1])\n", "targets_options['Zx_upper'] = np.array([0.6, 1.1])\n", "targets_options['Rs_lower_outer'] = np.array([0.9, -1.95])\n", "targets_options['Rs_upper_outer'] = np.array([0.9, 1.95])\n", "\n", "\n", "_, target_values = VCs.calculate_targets(eq, targets, targets_options)\n", "\n", "# print\n", "for i in range(len(targets)):\n", " print(targets[i] + \" = \" + str(target_values[i]))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also define **custom** shape targets by following the template in the next cell.\n", "\n", "For example, here we show how to calculate the $R_{gap}$ distance in the lower divertor. This is defined (in MAST-U at least) as the radial position of the point on the separatrix at $Z=-1.5$ and the wall (on the right side). " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# any new target function should take the equilibrium object as input\n", "import shapely as sh # requires the shapely package\n", "\n", "def R_gap_lower(eq):\n", " \n", " # find contour object for psi_boundary\n", " if eq._profiles.flag_limiter:\n", " cs = plt.contour(\n", " eq.R, eq.Z, eq.psi(), levels=[eq._profiles.psi_bndry]\n", " )\n", " else:\n", " cs = plt.contour(\n", " eq.R, eq.Z, eq.psi(), levels=[eq._profiles.xpt[0][2]]\n", " )\n", " plt.close() # this isn't the most elegant but we don't need the plot iteq\n", "\n", " # for each item in the contour object there's a list of points in (r,z) (i.e. a line)\n", " psi_boundary_lines = []\n", " for i, item in enumerate(cs.allsegs[0]):\n", " psi_boundary_lines.append(item)\n", "\n", " # use the shapely package to find where each psi_boundary_line intersects the z line\n", " gaps = []\n", " z_line = np.array([[0.5,-1.5],[0.8,-1.5]])\n", " curve1 = sh.LineString(z_line)\n", " for j, line in enumerate(psi_boundary_lines):\n", " curve2 = sh.LineString(line)\n", "\n", " # find the intersection point(s)\n", " intersection = curve2.intersection(curve1)\n", "\n", " # extract intersection point(s)\n", " if intersection.geom_type == \"Point\":\n", " gaps.append(np.squeeze(np.array(intersection.xy).T))\n", " elif intersection.geom_type == \"MultiPoint\":\n", " gaps.append(\n", " np.squeeze(\n", " np.array([geom.xy for geom in intersection.geoms])\n", " )\n", " )\n", "\n", " gap_point = np.array(gaps).squeeze()\n", "\n", " return gap_point[0] # select R position\n", "\n", "# we then create a list of lists to store the target name and its function, e.g. [[\"new_target_name\", ...], [new_target_function(eq),...]]\n", "non_standard_targets = [[\"R_gap_lower\"], [R_gap_lower]]\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# we can now include our custom target functions\n", "all_target_names, all_target_values = VCs.calculate_targets(eq, targets, targets_options, non_standard_targets)\n", "\n", "# print\n", "for i in range(len(all_target_names)):\n", " print(all_target_names[i] + \" = \" + str(all_target_values[i]))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us visualise the targets on our equilibrium. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot the resulting equilbria \n", "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)\n", "ax1.grid(True, which='both')\n", "ax1.contour(eq.R, eq.Z, eq.psi(), levels=[eq.psi_bndry], colors='r')\n", "eq.tokamak.plot(axis=ax1, show=False)\n", "ax1.plot(eq.tokamak.wall.R, eq.tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "\n", "ax1.scatter(all_target_values[0], 0.0, s=100, color='green', marker='x', zorder=20, label=all_target_names[0])\n", "ax1.scatter(all_target_values[1], 0.0, s=100, color='blue', marker='x', zorder=20, label=all_target_names[1])\n", "ax1.scatter(all_target_values[2], all_target_values[3], s=100, color='m', marker='*', zorder=20, label=f\"({all_target_names[2]},{all_target_names[3]})\")\n", "ax1.scatter(all_target_values[4], all_target_values[5], s=100, color='k', marker='*', zorder=20, label=f\"({all_target_names[4]},{all_target_names[5]})\")\n", "ax1.scatter(all_target_values[6], -1.9 , s=100, color='k', marker='x', zorder=20, label=f\"{all_target_names[6]}\")\n", "ax1.scatter(all_target_values[7], 1.95 , s=100, color='r', marker='x', zorder=20, label=f\"{all_target_names[7]}\")\n", "ax1.scatter(all_target_values[8], -1.95 , s=100, color='orange', marker='o', zorder=5, label=f\"{all_target_names[8]}\")\n", "\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "ax1.legend(loc=\"center\")\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Calculating the VCs\n", "\n", "Now we've defined the targets we're interested in, we can begin calculating the **shape** and **virtual circuit** matrices. The following is a brief outline of how they're calculated:\n", "\n", "##### 1. Initialise and solve the base equilibrium\n", "Before computing derivatives, the static forward GS solver is run with some initial coil currents. Following this we store the:\n", "- equilibrium state `eq` and plasma profiles `profiles`.\n", "- plasma current vector $I_y$ (which defines the amount of plasma current at each computational grid point, restricted to the limiter region to save computational resources).\n", "- target values $\\vec{T} = [T_1,\\ldots,T_{N_T}]$ are evaluated.\n", "\n", "This establishes a reference state before we perturb the coil currents - we have already carried out these steps above (the $I_y$ vector is stored within `profiles`).\n", "\n", "##### 2. Find the appropriate coil current perturbations\n", "Next, we need to identify the appropriate coil current perturbations $\\vec{\\delta I}$ to use in the finite differences. To do this, we aim to scale the starting guess `starting_dI` (denoted $\\vec{\\delta I}^{start}$) according to the target tolerance on the plasma current vector (`target_dIy` = $\\delta I_y^{target}$). While the user can choose `starting_dI`, the default is given by\n", "\n", "$$ \\vec{\\delta I}^{start} := | \\vec{I} | \\times \\delta I_y^{target}. $$\n", "\n", "For each coil current $j \\in \\{1,\\ldots,N_c \\}$, this starting guess is scaled as follows:\n", "1. Perturb coil current $j$: \n", "\n", "$$ I_j^{new} := I_j + \\delta I_j^{start}.$$\n", "\n", "2. Solve the equilibrium with the updated $j$ th coil current (others are left unchanged) and store the plasma current vector $\\vec{I}_y^{new}$.\n", "\n", "3. The starting guess for the $j$ th coil current perturbation is then scaled by the relative norm of the plasma current change (and the predefined target tolerance `target_dIy`):\n", "\n", "$$ \\delta I_j^{final} = \\delta I_j^{start} \\times \\left( \\delta I_y^{target} \\frac{\\| \\vec{I}_y^{start} \\|}{\\| \\vec{I}_y^{new} - \\vec{I}_y^{start} \\|} \\right).$$\n", "\n", "If this relative norm is larger than $\\delta I_y^{target}$, then the starting guess $\\delta I_j^{start}$ needs to be made smaller (and vice versa). \n", "\n", "After this, we have our scaled perturbations $\\vec{\\delta I}^{final}$ ready to use in the finite differences. \n", "\n", "##### 3. Find the finite differences (i.e. the shape matrix)\n", "For each coil current $j \\in \\{1,\\ldots,N_c \\}$:\n", "\n", "1. Perturb coil current $j$: \n", "\n", "$$ I_j^{new} := I_j + \\delta I_j^{final}.$$\n", "\n", "2. Solve the equilibrium with the updated $j$ th coil current (others are left unchanged) and store the plasma current vector $\\vec{I}_y^{new}$.\n", "\n", "3. Using the new target values $\\vec{T}^{new}$ from the equilibrium, calculate the $j$ th column of the shape matrix:\n", "\n", "$$ S_{:,j} = \\frac{\\vec{T}^{new} - \\vec{T}}{\\delta I_j^{final}}. $$\n", "\n", "This gives the sensitivity of each of the targets $T_i$ with respect to a change in the $j$ th coil current.\n", "\n", "Note that we can also obtain the Jacobian matrix of the plasma current vector (using $\\vec{I}_y^{new}$) with respect to the coil currents: $\\frac{\\partial \\vec{I}_y}{\\partial \\vec{I}}$.\n", "\n", "##### 4. Find the virtual circuit matrix\n", "Once the full shape matrix $S \\in \\Reals^{N_T \\times N_c} $ is known, the **virtual circuit matrix** is computed as:\n", "\n", "$$ V = (S^T S)^{-1} S^T \\in \\Reals^{N_c \\times N_T}.$$\n", "\n", "This matrix provides a mapping from requested shifts in the targets to the shifts in the coil currents required." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# define which coils we wish to calculate the shape (finite difference) derivatives (and therefore VCs) for\n", "coils = eq.tokamak.coils_list[0:12]\n", "print(coils)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# here we'll look at a subset of the coils and the targets\n", "coils = ['D1', 'D2', 'D3', 'Dp', 'D5', 'D6', 'D7', 'P4', 'P5']\n", "targets = ['R_in', 'R_out', 'Rx_lower', 'Zx_lower', 'Rs_lower_outer']\n", "non_standard_targets = [[\"R_gap_lower\"], [R_gap_lower]]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# calculate the shape and VC matrices as follows\n", "VCs.calculate_VC(\n", " eq=eq,\n", " profiles=profiles,\n", " coils=coils,\n", " targets=targets,\n", " targets_options=targets_options,\n", " target_dIy=1e-3,\n", " non_standard_targets=non_standard_targets,\n", " starting_dI=None,\n", " min_starting_dI=50,\n", " verbose=True,\n", " VC_name=\"VC_for_lower_targets\", # name for the VC\n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# these are the finite difference derivatives of the targets wrt the coil currents\n", "shape_matrix = 1.0*VCs.VC_for_lower_targets.shape_matrix\n", "print(f\"Dimension of the shape matrix is {shape_matrix.shape} [no. of targets x no. of coils]. Has units [m/A]\")\n", "\n", "# these are the VCs corresponding to the above shape matrix\n", "VCs_matrix = 1.0*VCs.VC_for_lower_targets.VCs_matrix\n", "print(f\"Dimension of the virtual circuit matrix is {VCs_matrix.shape} [no. of coils x no. of targets]. Has units [A/m].\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### How do we make use of the VCs?\n", "\n", "Now that we have the VCs, we can use them to idenitfy the coil current shifts required to change the targets by a certain amount. \n", "\n", "For example, we will ask for shifts in a few shape targets and observe what happens to the equilibrium once we apply the new currents from the VCs. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# let's remind ourselves of the targets\n", "print(targets)\n", "for i in range(len(non_standard_targets[0])):\n", " print(non_standard_targets[0][i])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# let's set the requested shifts (units are metres) for the full set of targets\n", "all_requested_target_shifts = [0.01, -0.01, 0.0, 0.0, 0.0, 0.0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We specify a perturbation in the targets $\\vec{\\Delta T}$ (the shifts above) and calculate the change in coil currents required to achieve this by calculating:\n", "\n", "$$ \\vec{\\Delta I} = V \\vec{\\Delta T}. $$\n", "\n", "Using $\\vec{\\Delta I}$ we then perturb the coil currents in our original equilibrium, re-solve the static forward Grad-Shafranov (GS) problem, and return. This is all done in the following cell. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# we apply the VCs to the desired shifts, using the VC_object we just created\n", "eq_new, profiles_new, all_target_names, new_target_values, old_target_values = VCs.apply_VC(\n", " eq=eq,\n", " profiles=profiles,\n", " VC_object=VCs.VC_for_lower_targets,\n", " all_requested_target_shifts=all_requested_target_shifts,\n", " verbose=True,\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot the resulting equilbria \n", "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)\n", "ax1.grid(True, which='both')\n", "ax1.plot(eq.tokamak.wall.R, eq.tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "eq.tokamak.plot(axis=ax1, show=False)\n", "\n", "ax1.contour(eq.R, eq.Z, eq.psi(), levels=[eq.psi_bndry], colors='r')\n", "ax1.contour(eq.R, eq.Z, eq_new.psi(), levels=[eq_new.psi_bndry], colors='b', linestyles=\"--\")\n", "\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can then measure the accuracy of the VCs by observing the difference between the requested change in shape targets vs. those enacted by the VCs. \n", " \n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for i in range(len(all_target_names)):\n", " print(f\"Difference in {all_target_names[i]} = {np.round(new_target_values[i] - old_target_values[i],3)} vs. requested = {(all_requested_target_shifts)[i]}.\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the plots below, we can see the actual shifts by the VCs are almost exactly as those requested (for the targets with actual shifts). Note how the unshifted targets also shift under the VCs due to the nonlinear coupling between targets and coil currents. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plots\n", "rel_diff = np.abs(np.array(new_target_values) - np.array(old_target_values))/np.abs(old_target_values)\n", "\n", "fig1, ((ax1), (ax2)) = plt.subplots(2, 1, figsize=(12,12), dpi=80)\n", "\n", "ax1.grid(True, which='both', alpha=0.5)\n", "ax1.scatter(all_target_names, all_requested_target_shifts, color='red', marker='o', s=150, label=\"Requested\")\n", "ax1.scatter(all_target_names, np.array(new_target_values) - np.array(old_target_values), color='royalblue', marker='o', s=75, label=\"Actual\")\n", "ax1.set_xlabel(\"Target\")\n", "ax1.set_ylabel(\"Shift [m]\")\n", "ax1.legend()\n", "# ax1.set_ylim([-max(targets_shift+non_standard_targets_shift)*1.1, max(targets_shift+non_standard_targets_shift)*1.1])\n", "\n", "ax2.grid(True, which='both', alpha=0.5)\n", "ax2.scatter(all_target_names, rel_diff, color='red', marker='o', s=150, edgecolors='black', label=\"Requested\")\n", "ax2.set_xlabel(\"Target\")\n", "ax2.set_ylabel(\"Relative shift\")\n", "labels = ax2.get_xticklabels()\n", "ax2.set_yscale('log')\n", "ax2.set_ylim([1e-6, 1e-0])\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Alternative VCs\n", "\n", "Note that when the `VC_calculate` method is called, a `VC_name` can be provided. This will store the key data in a subclass with `VC_name`. This useful for recalling from which targets and coils a given shape matrix or VC matrix was calculated. \n", "\n", "Note that if no `VC_name` is provided the default is used (\"latest_VC\")." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# display the attributes\n", "print(VCs.VC_for_lower_targets.__dict__.keys())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can then follow up by calculating alternative VCs for different target and coil combinations. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# choose new coils and targets\n", "coils = ['D1', 'D2', 'D3', 'Dp', 'D5', 'D6', 'D7', 'P4', 'P5']\n", "targets = ['Rx_upper', 'Zx_upper', 'Rs_upper_outer'] # this time we use the upper targets\n", "targets_options = dict()\n", "targets_options['Rx_upper'] = np.array([0.6, 1.1])\n", "targets_options['Zx_upper'] = np.array([0.6, 1.1])\n", "targets_options['Rs_upper_outer'] = np.array([1.15, 2.1])\n", "\n", "\n", "# calculate the shape and VC matrices\n", "VCs.calculate_VC(\n", " eq=eq,\n", " profiles=profiles,\n", " coils=coils,\n", " targets=targets,\n", " targets_options=targets_options,\n", " target_dIy=1e-3,\n", " non_standard_targets=None,\n", " starting_dI=None,\n", " min_starting_dI=50,\n", " verbose=True,\n", " VC_name=\"VC_for_upper_targets\",\n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# display the attributes\n", "print(VCs.VC_for_upper_targets.__dict__.keys())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now apply this VC using the `apply_VC` method." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# we'll shift the upper strikepoint (R value) \n", "all_requested_target_shifts = [0.0, 0.0, 0.02]\n", "\n", "# we apply the VCs to the some desired shifts, using the VC_object we just created\n", "eq_new, profiles_new, target_names, new_target_values, old_target_values = VCs.apply_VC(\n", " eq=eq,\n", " profiles=profiles,\n", " VC_object=VCs.VC_for_upper_targets,\n", " all_requested_target_shifts=all_requested_target_shifts,\n", " verbose=True,\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "VCs.VC_for_upper_targets.targets" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot the resulting equilbria \n", "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)\n", "ax1.grid(True, which='both')\n", "ax1.plot(eq.tokamak.wall.R, eq.tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "eq.tokamak.plot(axis=ax1, show=False)\n", "\n", "ax1.contour(eq.R, eq.Z, eq.psi(), levels=[eq.psi_bndry], colors='r')\n", "ax1.contour(eq.R, eq.Z, eq_new.psi(), levels=[eq_new.psi_bndry], colors='b', linestyles=\"--\")\n", "\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plots\n", "rel_diff = np.abs(np.array(new_target_values) - np.array(old_target_values))/np.abs(old_target_values)\n", "\n", "fig1, ((ax1), (ax2)) = plt.subplots(2, 1, figsize=(12,12), dpi=80)\n", "\n", "ax1.grid(True, which='both', alpha=0.5)\n", "ax1.scatter(target_names, all_requested_target_shifts, color='red', marker='o', s=150, label=\"Requested\")\n", "ax1.scatter(target_names, np.array(new_target_values) - np.array(old_target_values), color='royalblue', marker='o', s=75, label=\"Actual\")\n", "ax1.set_xlabel(\"Target\")\n", "ax1.set_ylabel(\"Shift [m]\")\n", "ax1.legend()\n", "# ax1.set_ylim([-max(targets_shift+non_standard_targets_shift)*1.1, max(targets_shift+non_standard_targets_shift)*1.1])\n", "\n", "ax2.grid(True, which='both', alpha=0.5)\n", "ax2.scatter(target_names, rel_diff, color='red', marker='o', s=150, edgecolors='black', label=\"Requested\")\n", "ax2.set_xlabel(\"Target\")\n", "ax2.set_ylabel(\"Relative shift\")\n", "labels = ax2.get_xticklabels()\n", "ax2.set_yscale('log')\n", "ax2.set_ylim([1e-6, 1e-0])\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "freegsnke4e", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.16" } }, "nbformat": 4, "nbformat_minor": 4 }